pinks = {'#860454','#c51b7c','#dc75ab','#f0b7da','#FF00FF','#ff81c0'};
greens = {'#276418','#4d921e','#7fbc42','#b9e084','#15b01a','#77AC30'};
purples = {'#410149','#762a84','#9b6fac','#c1a5cd','#800080','#7e1e9c'};
browns = {'#523107','#7f3c0a','#bf812c','#e2c17e','#a0522d','#A2142F'};
pools = {'#003d2e','#01665e','#379692','#81cdc1','#13eac9','#04d8b2'};
reds = {'#68011d','#b5172f','#d75f4e','#f7a580','#ff0000','#e50000'};
blues = {'#062e61','#2265ad','#4295c1','#93c5dc','#40e0d0','#00FFFF'};


figure(96)
clf
hold on
%% Linear polarization
Files_linear = {'data_20230330_100VR_30G_10kHz_REMPI_scan_21', ... %P2
    'data_20230330_100VR_30G_10kHz_REMPI_scan_22', ... % Q3
    'data_20230419_100VR_30G_10kHz_REMPI_scan_2', ... % Q3
    'data_20230328_100VR_30G_10kHz_REMPI_scan_19', ... % R3, this one is  actually polarization without waveplates 
    'data_20230330_100VR_30G_10kHz_REMPI_scan_24',... % Q2
    'data_20230331_100VR_30G_10kHz_REMPI_scan_2', ... % Q1
    'data_20230331_100VR_30G_10kHz_REMPI_scan_4', ... % R0
    'data_20230331_100VR_30G_10kHz_REMPI_scan_5', ... % R1
    'data_20230331_100VR_30G_10kHz_REMPI_scan_7'; %, ... % R2
    };

subplot(2,1,1)
hold on
for i = 1:length(Files_linear)
    d = load(Files_linear{i});
    x = d.x;
    ys = d.ys;
    eys = d.es;
    yb = d.yb;
    eyb = d.eb;
    errorbar(x,abs(ys-yb),sqrt(eys.^2+eyb.^2),'o','markersize',8,'linewidth',2,'color', hex2RGB(greens{1}),'markeredgecolor',hex2RGB(greens{1}),'markerfacecolor',hex2RGB(greens{4}))

    if i<2 || i>4
        xft = linspace(min(x),max(x),100);
        c0 = [100,50,mean(x)];
        lb = [0,0.1,mean(x)-20];
        ub = [5000,800,mean(x)+20];
        [c,cerr,fitval] = fit_wrapper(@(b,x)(b(1)*b(2).*2./((4*(x-b(3)).^2+b(2).^2)) ),c0,x,abs(ys-yb),sqrt(eys.^2+eyb.^2),xft,lb,ub);
        plot(xft,fitval,'-','linewidth',1.5,'color', hex2RGB(greens{1}))
    end
    
end

ylim([0 120])
% xlim([-11000,0])
box on; %grid on
xlabel('666 nm laser frequency 449850 GHz + x (MHz)')
ylabel('KRb ion counts')
set(gca,'fontsize',16)
set(gcf,'color','white')


% plot Rb
Files = {'data_20shots_20230521_100VR_30G_10kHz_REMPI_scan_3' ...
    'data_20shots_20230522_100VR_30G_10kHz_REMPI_scan_2'};




d1 = load(Files{1});
x1 = d1.x;
y1 = d1.yb;
ey1 = d1.eb;

d2 = load(Files{2});
x2 = d2.x-15;
y2 = d2.ys;
ey2 = d2.es;

x = [x1; x2];
y = [y1; y2];
ey = [ey1; ey2];

G = findgroups(x);
x = splitapply(@mean, x, G);
y = splitapply(@mean, y, G);
ey = splitapply(@(x) sqrt(mean(x.^2,1)),ey, G);


subplot(2,1,2)
hold on
errorbar(x,y,ey,'o','markersize',8,'linewidth',2,'color', hex2RGB(blues{1}),'markeredgecolor',hex2RGB(blues{1}),'markerfacecolor',hex2RGB(blues{4}))


xft = linspace(min(x),max(x),1000);
c0 = [300,50,70,  150,20,150,  150,20,255, 40,40,310,  40,40,355 ];
lb = [300-100,10,100-60,  100,5,180-40,  100,5,280-40,  2,10,340-40, 2,0.1,385-40];
ub = [80000,500,100,  1000,100,160, 1000,100,260, 100,100,315, 400,100,355];
[c,cerr,fitval] = fit_wrapper(@(b,x)(b(1)*b(2).*2./((4*(x-b(3)).^2+b(2).^2))+b(4)*b(5).*2./((4*(x-b(6)).^2+b(5).^2))+b(7)*b(8).*2./((4*(x-b(9)).^2+b(8).^2)) ...
     +b(10)*b(11).*2./((4*(x-b(12)).^2+b(11).^2))  +b(13)*b(14).*2./((4*(x-b(15)).^2+b(14).^2))),c0,x,abs(y),ey,xft,lb,ub);
plot(xft,fitval,'-','linewidth',1.5,'color', hex2RGB(blues{1}))


box on
xlim([0 430])
ylim([0 80])
xlabel('420 nm laser frequency 713288.4 GHz + x (MHz)')
ylabel('Rb ion counts')
set(gca,'fontsize',16)
set(gcf,'color','white')



ep = 1e-4;
% use B = 31.2G

w = 0.5*[0.71,0.51,0.45]; % plus, minus, pi
w = 0.5*[0.35/29.3*35.4,0.51/77.03*63.15,0.2]; % plus, minus, pi
load('Rb_REMPI_LinePred_sig_plus_31.mat')
x_m1_plus = mF_m1(find(mF_m1_dip>ep));
x_p1_plus = mF_p1(find(mF_p1_dip>ep));
x_p0_plus = mF_p0(find(mF_p0_dip>ep));
y_m1_plus = w(1)*mF_m1_dip(find(mF_m1_dip>ep));
y_p1_plus = w(1)*mF_p1_dip(find(mF_p1_dip>ep));
y_p0_plus = w(1)*mF_p0_dip(find(mF_p0_dip>ep));



load('Rb_REMPI_LinePred_sig_minus_31.mat')
x_m1_minus = mF_m1(find(mF_m1_dip>ep));
x_p1_minus = mF_p1(find(mF_p1_dip>ep));
x_p0_minus = mF_p0(find(mF_p0_dip>ep));
y_m1_minus = w(2)*mF_m1_dip(find(mF_m1_dip>ep));
y_p1_minus = w(2)*mF_p1_dip(find(mF_p1_dip>ep));
y_p0_minus = w(2)*mF_p0_dip(find(mF_p0_dip>ep));


load('Rb_REMPI_LinePred_pi_31.mat')
x_m1_pi = mF_m1(find(mF_m1_dip>ep));
x_p1_pi = mF_p1(find(mF_p1_dip>ep));
x_p0_pi = mF_p0(find(mF_p0_dip>ep));
y_m1_pi = w(3)*mF_m1_dip(find(mF_m1_dip>ep));
y_p1_pi = w(3)*mF_p1_dip(find(mF_p1_dip>ep));
y_p0_pi = w(3)*mF_p0_dip(find(mF_p0_dip>ep));



% display(mf_m1_dip)
% offset = min(mf_p1(abs(mf_p1_dip)>0));
match = 65;
% mf_p1m = mf_p1 - offset + match;
% mf_m1m = mf_m1 - offset + match;
% mf_p0m = mf_p0 - offset+ match;


stem(x_p1_pi + match,10*500*y_p1_pi ,'marker', 's','color','red','linewidth',2,'displayname','|11\rangle')
stem(x_p1_plus + match,10*500*y_p1_plus ,'marker', 's','color','red','linewidth',2,'displayname','|11\rangle')
stem(x_p1_minus + match,10*500*y_p1_minus ,'marker','s','color', 'red','linewidth',2,'displayname','|11\rangle')

stem(x_m1_pi + match,10*500*y_m1_pi ,'marker','o', 'color','green','linewidth',2,'displayname','|1,-1\rangle')
stem(x_m1_plus + match,10*500*y_m1_plus ,'marker','o','color', 'green','linewidth',2,'displayname','|1,-1\rangle')
stem(x_m1_minus + match,10*500*y_m1_minus ,'marker','o','color' ,'green','linewidth',2,'displayname','|1,-1\rangle')

stem(x_p0_pi + match,10*500*y_p0_pi ,'marker','d', 'color','blue','linewidth',2,'displayname','|10\rangle')
stem(x_p0_plus + match,10*500*y_p0_plus , 'marker','d','color', 'blue','linewidth',2,'displayname','|10\rangle')
stem(x_p0_minus + match,10*500*y_p0_minus ,'marker','d','color', 'blue','linewidth',2,'displayname','|10\rangle')




legend









function color = hex2RGB(hex)
color = sscanf(hex(2:end),'%2x%2x%2x',[1 3])/255;
end


function [c,cerr,fit_val] = fit_wrapper(func_hand,c0,x,y,ye,xft,lb,ub)
    if nargin < 5
        ye = 0;
    end
    if nargin < 6
        xft = x;
    end
    opts = optimset('Display','off');
    % [c,~,R,~,~,~,J] = lsqcurvefit(@gaussian2d,c0,xy,im(:),[],[],opts);
    if ye == 0
        [c,R,J] = nlinfit(x,y,func_hand,c0,opts); % slightly faster
    else
%         [c,R,J] = nlinfit(x,y,func_hand,c0,opts,'weights',1./ye.^2);
         [c,~,R,~,~,~,J] = lsqcurvefit(func_hand,c0,x,y,lb,ub);
    end
    ci = nlparci(c, R,'jacobian',J,'alpha',.3173); % 1 sigma
    cerr = abs(ci(:,2)-ci(:,1))'/2;
    fit_val = func_hand(c,xft);
end
